close all;clear all;clc;
% Parameters
A = 1.2;
g = 2;
lambda = 0.2;
eta = 0.6;

% Functions
tauf = @(g, gprime, lambda, eta, A, k)...
    1/(g-gprime)*(log(A)-log((exp((g-gprime)*k*eta)-1)/((g-gprime)*eta)...
    + (1-k)*((g/lambda)*exp((g-lambda-gprime)*k*eta)...
    - (g-lambda)/lambda*exp((g-gprime)*k*eta))));

alphaf = @(g, gprime, eta, A, k, tau)...
    (A*exp(gprime*(tau+k*eta))...
    - exp(g*tau)*(exp(g*k*eta)-exp(gprime*k*eta))/((g-gprime)*eta))/(1-k);

Abar = @(g, gprime, lambda, eta)...
    (exp((g-gprime)*eta/(lambda*eta)*log(g/(g-lambda)))-1)/((g-gprime)*eta);

k0f = @(g, gprime, lambda, eta, A)...
    fsolve(@(x) (exp((g-gprime)*x*eta)-1)/((g-gprime)*eta)...
    + (1-x)*((g/lambda)*exp((g-lambda-gprime)*x*eta)...
    - (g-lambda)/lambda*exp((g-gprime)*x*eta)) - A, 1);

A0f = @(g, gprime, lambda, eta, k)...
    ((exp((g-gprime)*k*eta)-1)/((g-gprime)*eta)...
    + (1-k)*((g/lambda)*exp((g-lambda-gprime)*k*eta)...
    - (g-lambda)/lambda*exp((g-gprime)*k*eta)))...
    /(exp((g-gprime)*k*eta)-(g-gprime)*eta...
    *((exp((g-gprime)*k*eta)-1)/((g-gprime)*eta)...
    + (1-k)*((g/lambda)*exp((g-lambda-gprime)*k*eta)...
    - (g-lambda)/lambda*exp((g-gprime)*k*eta))));

g0f = @(g, lambda, eta) g - lambda - log(1/(g*(eta-1/lambda*log(g/(g-lambda)))))...
    /(1/lambda*log(g/(g-lambda)));

hf = @(g, lambda, eta, k) lambda*exp(lambda*k*eta)/(exp(lambda*k*eta)-1);

kcf = @(g, gprime, lambda, eta) fsolve(@(k) (1-k)*g*eta + g/lambda...
    - (g-lambda)/lambda*exp(lambda*k*eta) - exp(-(g-lambda-gprime)*k*eta), 0);

knc = @(g, gprime, eta, A)...
    1/((g-gprime)*eta)*log(A*(g-gprime)*eta+1);

gprimebar = fsolve(@(x) 1/((g-x)*eta)*log(A*(g-x)*eta+1)...
    - 1/(lambda*eta)*log(g/(g-lambda)), 0);

% Calculations
step = 1e-3;
gprime = [0.1, 0.3, 0.5, 0.7];
k = NaN([length(gprime)...
    round(1/((g-max(gprime))*eta)*log(A*(g-max(gprime))*eta+1)/step)+1]);

for i = 1:length(gprime)
    kend(i) = floor(1/((g-gprime(i))*eta)*log(A*(g-gprime(i))*eta+1)*round(1/step))+1;
    k(i,1:kend(i)) = 0:step:1/((g-gprime(i))*eta)*log(A*(g-gprime(i))*eta+1);
    for j = 1:length(k(i,:))
        if ~isnan(k(i,j))
            tau(i,j) = max(tauf(g, gprime(i), lambda, eta, A, k(i,j)),0);
            tauPlusKeta(i,j) = tau(i,j) + k(i,j)*eta;
        else
            tau(i,j) = NaN;
            tauPlusKeta(i,j) = NaN;
        end
    end
    kc(i) = kcf(g, gprime(i), lambda, eta);
    A0(i) = A0f(g, gprime(i), lambda, eta, kc(i));
end

% Points for figures
for i = 1:length(gprime)
    for j = 1:length(k(i,:))
        if (j==1)|(j==length(k(i,:)))|(mod(j,100)==0)
            k_figure(i, j) = k(i, j);
            tauPlusKeta_figure(i, j) = tauPlusKeta(i, j);
        else
            k_figure(i, j) = NaN;
            tauPlusKeta_figure(i, j) = NaN;
        end
    end
end

% Figure
h = figure;
set(h, 'defaultLegendInterpreter','latex',...
    'defaultAxesTickLabelInterpreter','latex');
a = axes;

plot(k(1,:), tauPlusKeta(1, :),'-.k', 'linewidth', 2);
hold on;
plot(k(2,:), tauPlusKeta(2, :),'-.k', 'linewidth', 2);
plot(k(3,:), tauPlusKeta(3, :),':k', 'linewidth', 2);
plot(k(4,:), tauPlusKeta(4, :),':k', 'linewidth', 2);

p1 = plot(k_figure(1,:), tauPlusKeta_figure(1,:),'-.k*',...
    'linewidth',2,'MarkerSize',8);
p2 = plot(k_figure(2,:), tauPlusKeta_figure(2,:),'-.ko',...
    'linewidth',2,'MarkerSize',8);
p3 = plot(k_figure(3,:), tauPlusKeta_figure(3,:),':ks',...
    'linewidth',2,'MarkerSize',8);
p4 = plot(k_figure(4,:), tauPlusKeta_figure(4,:),':kd',...
    'linewidth',2,'MarkerSize',8);

xl = xlim();
yl = ylim();
set(a, 'xlim', [0 1], 'ylim', [0 yl(2)], 'XTick',[0 0.5 1],...
    'XTicklabel',{'$0$' '$0.5$' '$1$'}, 'YTick',[0.1 0.3 0.5],...
    'YTicklabel',{'$0.1$' '$0.3$' '$0.5$'}, 'fontsize', 16);
yl = ylim();
aux0 = xl(1):0.02:xl(2);
auy0 = aux0*eta;
plot(aux0, auy0, 'k.','MarkerSize',4, 'HandleVisibility','off');

text(0.2, 0.1, 'slope $=\eta$', 'interpreter','latex', 'Fontsize', 16);

xlabel('$k$', 'fontsize', 20, 'interpreter','latex');
ylabel('$\tau^{*}+k\eta$', 'fontsize', 20, 'interpreter','latex');
lgd = legend([p1, p2, p3, p4], ...
    {'$g^{\prime}=0.1$', '$g^{\prime}=0.3$', '$g^{\prime}=0.5$',...
    '$g^{\prime}=0.7$'},...
    'location', 'southeast', 'fontsize', 16);
lgd.ItemTokenSize(1) = 24;

%saveas(h, 'tauplusketa_gprime.eps');


